#ifndef WARRAY_IMPL_FFT_RADIX2_H_INCLUDED
#define WARRAY_IMPL_FFT_RADIX2_H_INCLUDED

#include <complex>
#include <cmath>
#include "../warray.h"

namespace warp {
namespace detail {

//! \brief Cooley-Tukey Radix 2 FFT implementation of a 1D DFT
template<typename T>
class cooley_tukey_radix2_1d_fft_implementor
{
public:
	cooley_tukey_radix2_1d_fft_implementor(size_t n)
	{
		// We can only do powers of two
#ifdef BOUNDS_CHECK
		assert( (n&(nn-1)) == 0 );
#endif
	}
	
	//! \brief Compute the 1D dft *not-in-place*
	void fwd(const std::complex<T>* l, int lstride, std::complex<T>* r, int rstride)
	{
		
	}
	
	//! \brief Compute the 1D dft *in-place*
	void fwd(std::complex<T>* r, int rstride)
	{
		T *data = reinterpret_cast<T*>(array);
		
		unsigned long n,mmax,m,j,istep,i;
		T wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;

		/*
		1 1
		3 9 swap(8,2) swap(9,3)
		5 5
		7 13 swap(6,12) swap(7,13)
		9 9
		11 5
		13 13
		15 9
		
		0 1 2 3 4 5 6 7
		0 1 4 3 2 5 6 7
		0 1 4 6 2 5 3 7
		*/

		size_t n = m_n << 1;
		j=1;
		for (i=1;i<n;i+=2) 
		{
			if (j > i) 
			{
				std::swap(data[j-1],data[i-1]);
				std::swap(data[j  ],data[i  ]);
			}
			m=nn;
			while (m >= 2 && j > m) 
			{
				j -= m;
				m >>= 1;
			}
			j += m;
		}
		mmax=2;
		
		while (n > mmax) 
		{
			istep=mmax << 1;
			theta=isign*(6.28318530717959/mmax);
			wtemp=sin(0.5*theta);
			wpr = -2.0*wtemp*wtemp;
			wpi=sin(theta);
			wr=1;
			wi=0;
			for (m=1;m<mmax;m+=2) 
			{
				for (i=m;i<=n;i+=istep) 
				{
					j=i+mmax;
					tempr=wr*data[j-1]-wi*data[j];
					tempi=wr*data[j]+wi*data[j-1];
					data[j-1]=data[i-1]-tempr;
					data[j]=data[i]-tempi;
					data[i-1] += tempr;
					data[i] += tempi;
				}
				wr=(wtemp=wr)*wpr-wi*wpi+wr;
				wi=wi*wpr+wtemp*wpi+wi;
			}
			mmax=istep;
		}
	}
};

} // namespace detail
} // namespace warp


#endif

